MD simulations and continuum theory of partially fluidized shear granular flows 
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We carry out a detailed comparison of soft particle molecular dynamics simulations with the 
theory of partially fluidized shear granular flows. We verify by direct simulations a constitutive 
relation based on the separation of the shear stress tensor into a fluid part proportional to the strain 
rate tensor, and a remaining solid part. The ratio of these two components is determined by the 
order parameter. Based on results of the simulations we construct the "free energy" function for the 
order parameter. We also present the simulations of the stationary deep 2D granular flows driven 
by an upper wall and compare it with the continuum theory. 
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When the ratio of shear to normal stress in a packed 
granular matter exceeds a certain threshold value, the 
granular matter yields and a flow ensues. In the last few 
years there have been many experimental, numerical and 
theoretical studies |, g, § |, g g g, g, |, [l(| 0, |l| that 
explored a broad range of granular flow conditions. While 
dilute granular flows can be well described by the kinetic 
theory of dissipative granular gases p3[ , dense granular 
flows still present significant difficulties for theoretical de- 
scription. A continuum theory of slow dense granular 
flows based on the so-called associated flow rule that re- 
lates the strain rate and the shear stress was proposed 
in Ref. This description neglects the effects of the 

dry friction between the grains and works only in a flu- 
idized state, so it cannot not describe hysteretic nature of 
the granular flow relevant for stick-slips avalanching 
0, etc. A similar model based on a Newtonian stress- 
strain constitutive relation with density dependent vis- 
cosity was proposed in Refs. ||. In this model, the 
viscosity diverges at the fluidization threshold when the 
density approaches the random close packing density of 
grains. 

Recently we proposed a phcnomcnological order pa- 
rameter (OP) description of the fluidization transition 
fj"5f . The OP specifies the ratio between solid and fluid 
parts of the stress tensor. The viscosity is defined as 
a ratio of the fluid part of the shear stress to the strain 
rate and remains finite at the fluidization threshold. This 
model yielded a good qualitative description of broad va- 
riety of phenomena occurring in granular flows. 

In this Letter we report on 2D soft particle molecular 
dynamics simulations performed to validate and quantify 
our OP theory. To fit the equation for the OP and stress- 
strain relation we performed simulations of the granular 
flow in a thin Couette geometry. The obtained set of 
equations was used to calculate the stress and velocity 
profiles in a different system, a thick granular layer under 
non-zero gravity driven by a moving heavy upper plate. 

The theory |l5| ] is based on a standard momentum 
conservation equation and the incompressibility condi- 



tion applicable for slow dense flows. To close the system, 
we assumed that the stress tensor a is comprised of two 
parts, a solid part a s , and a fluid part (taken in a 
purely Newtonian form) 



(1) 



where pf is the isotropic "partial" fluid pressure, pf is 
the viscosity coefficient associated with the fluid stress 
tensor. We set the fluid part of off-diagonal components 
of the stress tensor to be proportional to the off-diagonal 
components of the full stress tensor with the proportion- 
ality coefficient being a function of the OP p, 



°lx = q(p)<ry*; <Ty X = (i - q(p))v y 



(2) 



We choose a fixed range for the OP such that it is zero in 
a completely fluidized state and one in a completely static 
regime. Thus, the function q(p) has the property q(0) = 
1, q(l) = 0. In Refs.Jl5| for simplicity we postulated 
q(p) = 1 — p. A similar relationship can be introduced 
for the diagonal terms of the stress tensor 
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The dynamics of the OP was assumed to be relax- 
ational in nature and controlled by the generic Ginzburg- 
Landau equation, 



^ = DV*p 
Dt H 



F( P ,S) 



(5) 



Here D/Dt is the material derivative, D is the diffusion 
coefficient, F(p, S) is the derivative of the free energy 
density which has a quartic polynomial form to account 
for the bistability near the solid-fluid transition. Control 
parameter 5 was taken to be a linear function of (j) = 
max \ <J mn / <r nn \ . 

The OP p which plays a pivotal role in the theory, 
should be associated with the "microscopic" properties 
of the granular assembly. At any moment of time all con- 
tacts among the grains can be classified as either fluid-like 
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or solid-like. A contact is considered fluid-like if two par- 
ticles slide past each other or briefly collide, and is solid- 
like if two particles are jammed together for longer than 
a characteristic collision time. Here we postulate that 
the OP can be introduced as a ratio between space-time 
averaged numbers of solid contacts (Z s ) and all contacts 
(Z) within a sampling area, 

p(y) = {Zj/-(Zj. (6) 

where (£) and £ stand for averaging of £ in space and 
time respectively. This definition satisfies both limiting 
cases: when a granulate is in a solid state all contacts 
are stuck and p = 1; when it is strongly agitated (Z s ) is 
zero and (Z) is small but finite, therefore p = 0. Our OP 
is expected to be sensitive to the degree of fluidization. 
A small rearrangement of the force network may result 
in strong fluctuations of p, while the solid fraction and 
granular temperature remain virtually constant. This 
quantity is difficult to measure in experiments, however 
it can be found in soft-particle molecular dynamics. 

Molecular Dynamics Simulations The grains are 
assumed to be non-cohesive, dry, inelastic disk-like par- 
ticles. Two grains interact via normal and shear forces 
whenever they overlap. For the normal impact we em- 
ployed spring- dashpot model This model accounts 
for repulsion and dissipation; the repulsive component is 
proportional to the degree of the overlap and the velocity 
dependent damping component simulates the dissipation. 
The model for shear force is based upon technique devel- 
oped in pTsj ] . The motion of a grain is obtained by inte- 
grating the Newton's equations with forces and torques 
produced by interactions with neighbors and walls. For 
detailed discussion of the advantages and limitations of 
the employed model see Refs. || [lO], |ll). The compu- 
tational domain spans L x x L y area, and is periodic in 
horizontal direction x. The grain diameters are uniformly 
distributed around mean with relative width A r to avoid 
crystallization effects ||. The material parameters of 
grains were chosen similar to Ref. [[ll] . All quantities are 
normalized by an appropriate combination of the average 
particle diameter d, mass m, and gravity g. 

We studied a thin (50 x 10) granular layer sandwiched 
between two "rough plates" under fixed external pressure 
P ext and zero gravity conditions. The rough plates were 
simulated by two straight chains of large grains "glued" 
together. Two opposite forces Fi = — F2 were applied to 
the plates along the horizontal x axis to induce shear 
stress in the bulk. We started with weak forces and 
slowly ramped them up in small increments. After that 
we ramped the shear forces down until the granular layer 
was jammed again. At every "stop" we measured all 
stress components, strain rate, and the OP by averaging 
over the whole layer and over time. These simulations 
were carried out at several values of the external pres- 
sure P ex t- Figure |l| shows p as a function of the nor- 
malized shear stress 8 — —o~ vx / P ex t- With this normal- 
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FIG. 1: p vs 5 for P ext = 20 (D),30 (O) and 40 (A). Open 
symbols correspond to ramp-up, and filled symbols to ramp- 
down in shear stress. Dashed line shows fit Eq. ^. Inset: 
density v vs p for the same three Pept- 
ization, the results of several simulations with different 
Pext collapse on a single bifurcation curve. We observe 
a quiescent state p = 1 until 8 reaches a certain critical 
value 81 » 0.32. This value differs slightly for differ- 
ent runs because of the finite system size and absence of 
self-averaging in the static regime. Above Si, p abruptly 
drops to approximately 0.15. At larger 5, the OP rapidly 
approaches zero. The return curve corresponding to the 
diminishing of the shear stress follows roughly the same 
path, and then continues to another (smaller) value of 
the shear stress, 82 ~ 0.23. At this point the OP jumps 
back to one, and the granular layer returns to a jammed 
state. The striking feature of this bifurcation diagram 
is the hysteretic behavior of the OP as a function of the 
shear stress. This hysteresis was anticipated in our model 
|]l5f , and now we are in a position to describe it quantita- 
tively. Assuming that there is an (unobserved) unstable 
branch of the bifurcation curve which merges with the 
stable branch at i5 « Si, we make a simple analytic fit, 

F(p, 8) = (l-p) (p 2 - 2p*p + pi exp[-A(S 2 Si)}) (7) 

with = 0.6, A = 25, 8* = 0.26 (see Figure |], dashed 
line) and use it in equation (0). The inset of Figure 1 
depicts the density v vs. the OP p for the same runs. 
The density stays almost constant in a wide range of the 
OP 0.1 < p < 1. This shows that unlike the particle 
density, our OP is a sensitive characteristic of slow dense 
granular flows reflecting subtle changes in the contact 
network and the structure of the stress distribution. 

We also probed the relaxation dynamics of the OP by 
studying the response of the system on small perturba- 
tion in the the hysteretic region, 5% < S < Si From 
these simulations we find that the intrinsic time scale of 
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FIG. 2: Ratios of the fluid stress components to the full 
stress components Oy/cry vs. p for P e;c t = 20 (□),30 (O) 
and 40 (A), a - shear stress components a yx (open symbols), 
o yx (filled symbols), line is a fit q(p) = (1 — p) 2 ' 5 , b - normal 
stress components a xx (closed symbols), a yy (open symbols), 
lines are the fits q x {p) = (1 - p) 1 ' 9 , 1y{p) = (1 — p 1 ' 2 ) 1,9 . 



the OP relaxation is rather small, 0(1). Our thin Cou- 
ette flow system did not allow us to probe the local cou- 
pling of the OP since p w const throughout the system. 
In the absence of such data this coupling was modeled 
by the linear diffusion term in (|5|) with D = const. 

The constitutive relation was fitted using the same 
Couette flow simulations. We analyzed the fluid stress 
o~fj and the solid stress of^ separately during our ramp- 
down simulations at three different values of P. Figure 
H shows the ratios of fluid tensor components to the cor- 
responding full tensor components as functions of p. We 
observe that data for <J s xy jo xy and <J yx jo yx for different P 
fall onto a single curve which is fitted by q(p) — (1 — p) 2 ' 5 
(Figure |[a). Repeating the procedure with diagonal el- 
ements yields different scaling, see Figure ||, b. A small 
but noticeable difference is evident between o~l x / o xx and 
a vv/ (J yy Detailed analysis shows that in fact fluid parts 
of the diagonal components of the stress tensor a' xx and 
crf y are nearly identical [fl7| , and the difference is mainly 
due to the solid part of the normal stresses. Functions 
q x , y (p) approach 1 as p — > 0, but they may have differ- 
ent functional form to reflect the anisotropy of the solid 
stress tensor. In our Couette flow, a xx and o~ yy can be 
fitted by q x (p) « (1 - p) 19 and q y (p) « (1 - p 12 ) 19 , 
respectively, see Figure ||,b. We observe that even in a 
partially fluidized regime, the "fluid phase" component 
behaves as a fluid with a well-defined isotropic "partial" 
pressure p / which is zero in a solid state and is becoming 
the full pressure in a completely fluidized state. 

To test the stress-strain relation (|l]) we plot — a^ yx vs 
^fyx, see Figure |3|. At small j yx all curves are close to 



FIG. 3: Stress-strain rate relation for a thin granular Couette 
flow at P ext = 20 (D),30 (O) and 40 (A). Fluid shear stress 
vs strain rate, the straight line is a constant viscosity fit = 
I27; Inset: full shear stress vs strain rate 



the same straight line a' yx « \2^ yx , i.e. the Newtonian 
scaling for fluid shear stress holds reasonably well. The 
deviations at large ^ yx are evidently caused by variations 
of density and temperature in the dilute regime as to 
be expected from kinetic theory of dilute granular flows 
ifLsf . The full shear stress, of course, does not vanish 
as j yx — ► (Figure |^, inset). Therefore, the standard 
viscosity coefficient defined as the ratio of the full shear 
stress and strain rate diverges at the fluidization thresh- 
old as observed in Ref . ||] . 

We applied our theoretical description which was for- 
mulated above on the basis of MD simulations of a thin 
Couette flow with no gravity, to a different system, a 
shear granular flow in a thick granular layer under grav- 
ity driven by the upper plate which was pulled with con- 
stant speed V (or constant force F). A similar system 
has been studied experimentally Refs. ^ ^|. We sim- 
ulated up to 20,000 particles in a periodic rectangular 
box under a heavy plate which was moved horizontally 
in x-direction. We systematically carried out comparison 
between MD simulations and the continuum theory us- 
ing the stress-strain relations and the specific form of the 
OP equation described above. We used no-flux boundary 
conditions for the OP and no-slip condition for the veloc- 
ity at the bottom plate. The boundary condition at the 
top plate where the flow can be in a dilute regime, is a 
separate issue which we do not address here. We assumed 
that the shear stress component in the bulk is specified 
by applied force, and calculated the velocity profile us- 
ing the constitutive relation (|l|)-(0). We found that the 
constitutive relations determined from the thin Couette 
experiment hold for this system as well. Selected results 
for the OP and velocity profiles are presented in Figure 
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mersed" in the solid phase. The fluid phase behaves as 
a simple Newtonian fluid for small shear rates when the 
density is almost constant, but exhibits shear thinning 
at larger shear rates when the density begins to drop. 
This regime can be described by the generalization of the 
theory Q which includes the equations for density and 
granular temperature following from the kinetic theory of 
dilute granular gases. Our simulations were limited to 2D 
geometry. While we anticipate that the structure of the 
model should remain unchanged in 3D systems, the spe- 
cific form of the fitting functions may vary. The authors 
are indebted to J. Gollub, P. Cvitanovic, T. Halsey, B. 
Beringer for useful discussions, and to J.C. Tsai for shar- 
ing his unpublished experimental data. This work was 
supported by the Office of the Basic Energy Sciences at 
the US DOE, grants W-31-109-ENG-38, and DE-FG03- 
95ER14516. Simulation were performed at the National 
Energy Research Scientific Computing Center. 



FIG. 4: Vertical profiles of p and Vx in a thick granular layer 
driven at the surface by a heavy plate for 5,000 particles, in 
the box L x = 50, L y = 100, a - P ext = 50, pulling velocity 
V = 5, D = 5, b - Pext = 10, pulling velocity V = 50, 
D = 1. Lines show the theoretical results obtained from the 
continuum model (|s]), ([?]), open symbols show numerical data. 

||. As seen from the Figure, the vertical profiles of the 
OP and the horizontal velocities are reasonably well de- 
scribed by the theory. However, for low pressure the hor- 
izontal velocity profiles deviate from the numerical data 
for low pressure runs, apparently because the viscosity 
coefficient is no longer a constant in a dilute region near 
the top plate. The only fitting parameter used was the 
diffusion constant D in the OP equation, which has not 
been determined in our simulations of the thin layer. It 
appears that the diffusion coefficient depends on the ap- 
plied pressure and the strain rate, however more detailed 
numerical experiments are needed. 

In conclusion, we calibrated the theory of partially 
granular flows |lq] on the basis of a series of 2D soft parti- 
cle molecular dynamics simulations. The OP which con- 
trols the fluidization transition, was defined as the frac- 
tion of solid- like contacts among particles. Measurements 
of the OP, the stress tensor, and the strain rate in a thin 
Couette cell allowed us to quantify the constitutive rela- 
tions based on the relaxational OP dynamics. We studied 
the flow structure of a thick surface driven granular flow 
under gravity and found the model predictions to be in a 
good quantitative agreement with soft-particle MD simu- 
lations. Our results support an intriguing interpretation 
for the OP description of dense and slow granular flows. 
The granular matter under shear stress appears to be 
similar to a multi-phase system with the fluid phase "im- 
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